/*
  rm15dump v 1.0 - Downloads data from Geoscan Research RM15-D resitivity meters.
  Copyright (C) 2010 John Donovan - mersey.viking@gmail.com

  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  
  http://www.apache.org/licenses/LICENSE-2.0
  
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/

/** \file main.cpp */

/**
\mainpage rm15dump

Downloads data from Geoscan Research (http://www.geoscan-research.co.uk/)
RM15-D resitivity meters, and outputs it into a variety of useful formats.
It has only been tested on non-multiplexed RM15-D units, but other hardware
(i.e. RM4 with DL10 and RM15) should also work. Any successes, failures,
or bugs in this software should be reported to the author.

\author  John Donovan
\date
2010
\version 1.0

Example usage:

./rm15dump -g"1,2" -w 30 -l 30 -p "/dev/ttyS0" --xyz --tiff outfile

Downloads the first two tiles of data from the RM15-D on serial port /dev/ttyS0,
sets the width and length of the tile to 30 metres, and outputs an XYZCSV file
(outfile.xyz) for importing into a third-party application, and a GeoTiff
(outfile.tif) file for importing into a GIS application. Use rm15dump --help
for information on all command line parameters.

\note
In the Geoscan documentation, what is referred to as a \a grid is called in
this code and documentation a \a tile. This is to distinguish the overall site
grid in which each tile has a place.

\note
Tile origin is always the lower left corner.

\note
To compile the code, you need Boost 1.37.0 (http://www.boost.org/),
Qt >= 4.2.0 (http://www.trolltech.com), 
argtable2 (http://argtable.sourceforge.net/), and GDAL 1.6.0
(http://www.gdal.org/) development libaries, and cmake 2.6
(http://www.cmake.org/).

\note
On Microsoft Windows machines, the serial port will be called COMn, where n is
a number usually from 1 to 4, but can be higher.

Licensed under the Apache License, Version 2.0
Copyright (C) 2010 John Donovan - mersey.coder@gmail.com
*/

#include <math.h>
#include <float.h>

#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>

#include <argtable2.h>

#include <gdal_priv.h>
#include <ogr_spatialref.h>

#include <QRegExp>
#include <QStringList>
#include <QTextStream>

#include "rm15d.h"

/**
\brief Specifies the coordinate of a tile of data in relation to the site's grid system.
*/
struct GridCoordinate {
  size_t _x; /*!< \brief Grid x coordinate of the tile. */
  size_t _y; /*!< \brief Grid y coordinate of the tile. */
};

enum Corners{TOP_LEFT, TOP_RIGHT, BOTTOM_LEFT, BOTTOM_RIGHT};

typedef std::pair<int, GridCoordinate> grid_coord_pair_t;
typedef std::map<int, GridCoordinate> grid_coord_map_t;
typedef grid_coord_map_t::iterator grid_coord_map_iter_t;

/** \brief Program options as they come from the command line. */
struct Options {
  std::vector<size_t> _tileIndexes; /*!< \brief Tile index as it comes from the RM15-D. Starts from 1. */
  grid_coord_map_t _gridPositions;  /*!< \brief For any tile that has been given a grid reference. */
  size_t _xsize;                    /*!< \brief Width, in metres, of a tile, i.e. the traverse length.*/
  size_t _ysize;                    /*!< \brief Length, in metres, of a tile. */
  size_t _xsamp;                    /*!< \brief Number of samples per traverse. */
  size_t _ysamp;                    /*!< \brief Number of traverses per metre. */
  std::string _prefix;              /*!< \brief Filename prefix. */
  size_t _baud;                     /*!< \brief Baud rate setting on the RM15-D. */
  std::string _port;                /*!< \brief Serial port name. */
  bool _tiff;                       /*!< \brief Should we generate a GeoTiff file? */
  bool _raw;                        /*!< \brief Should we dump the raw data from the RM15-D to a file? */
  bool _xyz;                        /*!< \brief Should we generate an XYZCSV file for importing into Geoplot? */
  bool _zigzag;                     /*!< \brief Traverse type, either parallel or zigzag. */
  size_t _frequency;                /*!< \brief Frequency setting of the RM15-D, in Hz. */
  std::string _comment;             /*!< \brief A comment to be put into the metadata. */
  double _nodata;                   /*!< \brief The no-data value to be used in a GeoTiff. */
  double _probeSeparation;          /*!< \brief The separation of the mobile probes for calculating resistivity. */
  double _originX;                  /*!< \brief The x-coordinate of the origin, in the chosen projection, of the tile. */
  double _originY;                  /*!< \brief The y-coordinate of the origin, in the chosen projection, of the tile. */
  double _axisX;                    /*!< \brief The x-coordinate of the x-axis, in the chosen projection, of the tile. */
  double _axisY;                    /*!< \brief The y-coordinate of the x-axis, in the chosen projection, of the tile. */
  std::string _coordinateRefSystem; /*!< \brief The projection in either EPSG code, proj.4 format, or WKT. */
  size_t _falseOriginX;             /*!< \brief The false origin offset to ensure all grid references are positive. Specified in tiles. */
  size_t _falseOriginY;             /*!< \brief The false origin offset to ensure all grid references are positive. Specified in tiles. */
  Corners _corner;                  /*!< \brief The corner from which the tile was started. */
  bool _clockwise;                  /*!< \brief The direction of the first traverse. */
};

/** \brief Global variable that stores all the command line options. */
Options opts;

/** \brief Stores the overall bounds of the data, in grid squares. */
struct GridBounds {
	GridBounds() : _gridXMin(INT_MAX), _gridXMax(0), _gridYMin(INT_MAX), _gridYMax(0) {}
	size_t _gridXMin;
	size_t _gridXMax;
	size_t _gridYMin;
	size_t _gridYMax;
};

GridBounds gridBounds;

const char* sourceXML = "<SimpleSource>\
<SourceFilename relativeToVRT=\"1\">%s</SourceFilename>\
<SourceBand>1</SourceBand>\
<SourceProperties RasterXSize=\"%d\" RasterYSize=\"%d\" DataType=\"Float32\" BlockXSize=\"%d\" BlockYSize=\"1\"/>\
<SrcRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\
<DstRect xOff=\"%d\" yOff=\"%d\" xSize=\"%d\" ySize=\"%d\"/>\
</SimpleSource>";

// Prototypes.
bool parseOptions(int argc, char* argv[]);
int parseGrids(struct arg_str* parent);
int parseTraverse(struct arg_str* parent);
int parseBaud(struct arg_int* parent);
void fixTiles(std::vector<Sample*>& sampleVec);
void mirrorTile(std::vector<Sample*>::iterator begin, std::vector<Sample*>::iterator end, size_t samplesPerTraverse);
void flipTile(std::vector<Sample*>::iterator begin, size_t samplesPerTraverse, size_t traversesPerTile);

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
int main(int argc, char* argv[]) {
  if (!parseOptions(argc, argv)) {
    return 1;
  }

  const size_t samplesPerTraverse = opts._xsize * opts._xsamp;
  const size_t samplesPerTile = samplesPerTraverse * opts._ysize * opts._ysamp;
  const size_t totalSamplesToRead = (*opts._tileIndexes.rbegin()) * samplesPerTile;
  
  std::vector<Sample*> sampleVec;
  RM15Device device(opts._port, opts._baud, opts._probeSeparation);
  int res = device.downloadData(sampleVec, totalSamplesToRead);
  
  if (res != 0) {
    return res;
  }
  
  fixTiles(sampleVec);
  
  GDALDriver* tiffDriver = 0;
  GDALDriver* vrtDriver = 0;
  GDALDataset* tiffDataSet = 0;
  GDALDataset* vrtDataSet = 0;

  if (opts._tiff) {
    GDALRegister_GTiff();
    GDALRegister_VRT();
    tiffDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
    vrtDriver = GetGDALDriverManager()->GetDriverByName("VRT");

    vrtDataSet = vrtDriver->Create((opts._prefix + ".vrt").c_str(), 
     ((gridBounds._gridXMax - gridBounds._gridXMin) + 1) * opts._xsize * opts._xsamp, 
     ((gridBounds._gridYMax - gridBounds._gridYMin) + 1) * opts._ysize * opts._ysamp, 
     1, GDT_Float32, 0);

    if (vrtDataSet) {
      vrtDataSet->GetRasterBand(1)->SetColorInterpretation(GCI_GrayIndex);
      vrtDataSet->GetRasterBand(1)->SetNoDataValue(opts._nodata);
    }
  }

  for (size_t tile = 0; tile < opts._tileIndexes.size(); ++tile) {
    char curTileStr[32];
    sprintf(curTileStr, "_%d", opts._tileIndexes[tile]);
    const size_t firstSample = (opts._tileIndexes[tile] - 1) * samplesPerTile;
    const size_t lastSample = opts._tileIndexes[tile] * samplesPerTile;

    // Write out metadata.
    {
      std::string fn = opts._prefix;
      fn += curTileStr;
      fn += ".log";
      FILE* f = fopen(fn.c_str(), "w");

      if (!f) {
        std::cerr << "ERROR: Unable to create file \"" << fn << "\"" << std::cerr;
        continue;
      }

      grid_coord_map_iter_t it = opts._gridPositions.find(opts._tileIndexes[tile]);

      if (it != opts._gridPositions.end()) {
        fprintf(f, "index %d %d\n", it->second._x, it->second._y);
      }

      fprintf(f, "tile width metres %d\n", opts._xsize);
      fprintf(f, "tile length metres %d\n", opts._ysize);
      fprintf(f, "samples per metre %d\n", opts._xsamp);
      fprintf(f, "traverses per metre %d\n", opts._ysamp);

      if (opts._zigzag) {
        fprintf(f, "traverse type zigzag\n");
      } else {
        fprintf(f, "traverse type parallel\n");
      }

      fprintf(f, "frequency %d\n", opts._frequency);
      fprintf(f, "nodata %f\n", opts._nodata);
      fprintf(f, "probe separation metres %f\n", opts._probeSeparation);
      fprintf(f, "comment\n%s\n", opts._comment.c_str());
      fclose(f);
    }

    if (opts._tiff) {
      std::string fn = opts._prefix;
      fn += curTileStr;
      fn += ".tif";

      char** papszOptions = NULL;
      char buf[32];
      sprintf(buf, "%d", opts._xsize * opts._xsamp);
      papszOptions = CSLSetNameValue(papszOptions, "BLOCKXSIZE", buf);
      papszOptions = CSLSetNameValue(papszOptions, "BLOCKYSIZE", "1");

      // It is usual to display high resistivity as black.
      papszOptions = CSLSetNameValue(papszOptions, "PHOTOMETRIC", "MINISWHITE");
      tiffDataSet = tiffDriver->Create(fn.c_str(), opts._xsize * opts._xsamp, opts._ysize * opts._ysamp, 1, GDT_Float32, papszOptions);

      if (tiffDataSet) {
        grid_coord_map_iter_t it = opts._gridPositions.find(opts._tileIndexes[tile]);
        double georef[6] = {0.0, 1.0 / opts._xsamp, 0.0, opts._ysize, 0.0, -1.0 / opts._ysamp};
        float* data = new float[opts._xsize * opts._xsamp * opts._ysize * opts._ysamp];

        if (it != opts._gridPositions.end()) {
          georef[0] = it->second._x * opts._xsize;
          georef[3] = (it->second._y + 1) * opts._ysize;
        }

        tiffDataSet->SetGeoTransform(georef);

        double min = DBL_MAX;
        double max = -DBL_MIN;
        double mean = 0.0;
        double stdDev = 0.0;

        for (size_t y = 0; y < opts._ysize * opts._ysamp; ++y) {
          for (size_t x = 0; x < samplesPerTraverse; ++x) {
            size_t index = firstSample;

            index += (y * samplesPerTraverse) + x;

            double res = 0.0;

            res = sampleVec[index]->value();

            if (res < min) {
              min = res;
            }

            if (res > max) {
              max = res;
            }

            mean += res / (opts._ysize * opts._ysamp * samplesPerTraverse);

            data[(y * samplesPerTraverse) + x] = static_cast<float>(res);
          }
        }

        for (size_t i = 0; i < opts._ysize * opts._ysamp * samplesPerTraverse; ++i) {
          stdDev += (mean - data[i]) * (mean - data[i]) / ((opts._ysize * opts._ysamp * samplesPerTraverse) - 1.0);
        }

        stdDev = sqrt(stdDev);
        char stddevStr[32];
        sprintf(stddevStr, "%f", stdDev);
        tiffDataSet->SetMetadataItem("STDDEV", stddevStr);

        GDALRasterBand* band = tiffDataSet->GetRasterBand(1);
        band->SetNoDataValue(opts._nodata);
        band->SetColorInterpretation(GCI_GrayIndex);
        band->RasterIO(GF_Write, 0, 0, opts._xsize * opts._xsamp, opts._ysize * opts._ysamp, data, opts._xsize * opts._xsamp, opts._ysize * opts._ysamp, GDT_Float32, 0, 0);

        if (vrtDataSet) {
          char buff[4096];
          sprintf(buff, sourceXML, fn.c_str(), tiffDataSet->GetRasterXSize(), tiffDataSet->GetRasterYSize(), tiffDataSet->GetRasterXSize(), 
            tiffDataSet->GetRasterXSize(), tiffDataSet->GetRasterYSize(), 
            (it->second._x - gridBounds._gridXMin) * opts._xsize * opts._xsamp, (gridBounds._gridYMax - it->second._y) * opts._ysize * opts._ysamp,
            tiffDataSet->GetRasterXSize(), tiffDataSet->GetRasterYSize());

          vrtDataSet->GetRasterBand(1)->SetMetadataItem("source_0", buff, "new_vrt_sources");
        }

        GDALClose(tiffDataSet);
        delete[] data;
      } else {
        std::cerr << "ERROR: Unable to create file \"" << fn << "\"" << std::cerr;
      }
    }

    if (opts._raw) {
      std::string fn = opts._prefix;
      fn += curTileStr;
      fn += ".raw";
      FILE* f = fopen(fn.c_str(), "w");

      if (!f) {
        std::cerr << "ERROR: Unable to create file \"" << fn << "\"" << std::cerr;
        continue;
      }

      for (size_t s = firstSample; s < lastSample; ++s) {
        fprintf(f, "%s\n", static_cast<const char*>(sampleVec[s]->raw().c_str()));
      }
      fclose(f);
    }

    if (opts._xyz) {
      std::string fn = opts._prefix;
      fn += curTileStr;
      fn += ".xyz";
      FILE* f = fopen(fn.c_str(), "w");

      if (!f) {
        std::cerr << "ERROR: Unable to create file \"" << fn << "\"" << std::cerr;
        continue;
      }

      for (size_t y = 0; y < opts._ysize * opts._ysamp; ++y) {
        for (size_t x = 0; x < samplesPerTraverse; ++x) {
          size_t index = firstSample;

          index += (y * samplesPerTraverse) + x;

          fprintf(f, "%d,%d,%0.6f\n", x, y, sampleVec[index]->value());
        }
      }
      fclose(f);
    }
  }

  if (vrtDataSet) {
    if (!(opts._coordinateRefSystem.empty() || opts._originX == 0.0 || opts._originY == 0.0 || opts._axisX == 0.0 || opts._axisY == 0.0)) {
      OGRSpatialReference srs;

      if (srs.SetFromUserInput(opts._coordinateRefSystem.c_str()) == OGRERR_NONE) {

        char* srsWKT = NULL;
        GDAL_GCP* pasGCPs = static_cast<GDAL_GCP*>(CPLMalloc(sizeof(GDAL_GCP) * 4));
        GDALInitGCPs(4, pasGCPs);

        const double x = opts._xsize * opts._xsamp;
        const double y = vrtDataSet->GetRasterYSize();

        // Transform the bottom left corner of the image.
        const double baseLen = sqrt(((opts._axisX - opts._originX) * (opts._axisX - opts._originX)) + ((opts._axisY - opts._originY) * (opts._axisY - opts._originY)));
        const double baseVecX = (opts._axisX - opts._originX) / baseLen;
        const double baseVecY = (opts._axisY - opts._originY) / baseLen;

        const double offsetX = (baseVecX * (static_cast<double>(gridBounds._gridXMin) - opts._falseOriginX) * opts._xsize) - (baseVecY * (static_cast<double>(gridBounds._gridYMin) - opts._falseOriginY) * opts._ysize);
        const double offsetY = (baseVecX * (static_cast<double>(gridBounds._gridYMin) - opts._falseOriginY) * opts._ysize) + (baseVecY * (static_cast<double>(gridBounds._gridXMin) - opts._falseOriginX) * opts._xsize);

        // Creates GCPs for the bottom left tile of the image.
        const double llx = opts._originX + offsetX;
        const double lly = opts._originY + offsetY;
        const double lrx = baseVecX * opts._xsize;
        const double lry = baseVecY * opts._ysize;
        const double ulx = -baseVecY * opts._xsize;
        const double uly = baseVecX * opts._ysize;
        const double urx = lrx + ulx;
        const double ury = lry + uly;

        // Top left.
        pasGCPs[0].dfGCPPixel = 0.0;
        pasGCPs[0].dfGCPLine = y - (opts._ysize * opts._ysamp);
        pasGCPs[0].dfGCPX = llx + ulx;
        pasGCPs[0].dfGCPY = lly + uly;

        // Top right.
        pasGCPs[1].dfGCPPixel = x;
        pasGCPs[1].dfGCPLine = y - (opts._ysize * opts._ysamp);
        pasGCPs[1].dfGCPX = llx + urx;
        pasGCPs[1].dfGCPY = lly + ury;

        // Bottom right.
        pasGCPs[2].dfGCPPixel = x;
        pasGCPs[2].dfGCPLine = y;
        pasGCPs[2].dfGCPX = llx + lrx;
        pasGCPs[2].dfGCPY = lly + lry;

        // Bottom left.
        pasGCPs[3].dfGCPPixel = 0.0;
        pasGCPs[3].dfGCPLine = y;
        pasGCPs[3].dfGCPX = llx;
        pasGCPs[3].dfGCPY = lly;

        srs.exportToWkt(&srsWKT);
        const char* srsWKTStart = srsWKT;
        vrtDataSet->SetProjection(srsWKTStart);
        // Uses the special case for four corners.
        double georef[6];
        GDALGCPsToGeoTransform(4, pasGCPs, georef, false);
        vrtDataSet->SetGeoTransform(georef);

        GDALDeinitGCPs(4, pasGCPs);
        CPLFree(pasGCPs);
        CPLFree(srsWKT);
      }
    }

    GDALClose(vrtDataSet);
  }
  
  // Clean up.
  for (size_t i = 0; i < sampleVec.size(); ++i) {
    delete sampleVec[i];
  }

  return EXIT_SUCCESS;
}

///////////////////////////////////////////////////////////////////////////////
/// The format is a string of comma-separated tile indexes optionally followed
/// by a +ve grid reference. The tile indexes can be specified in any order.
///
/// E.g. 1,2,3
///
/// 1[0:0],2[1:0]
///
/// 1,3[5:4],2
///////////////////////////////////////////////////////////////////////////////
int parseGrids(struct arg_str* parent) {
  QString grids(parent->sval[0]);
  grids.remove(' ');
  
  // Regular expression that matches tile references from the command line.
  QRegExp gridPattern("(\\d+)(?:,|$)");
  // Regular expression that matches a tile index and grid position.
  QRegExp coordPattern("(\\d+)\\[(\\d+):(\\d+)\\](?:,|$)");


  int pos = 0;
  while ((pos = gridPattern.indexIn(grids, pos)) != -1) {
    int gridIndex = gridPattern.cap(1).toInt();
    if (gridIndex < 1) {
      std::cerr << "WARNING: Ignoring tile with index of " << gridIndex << std::endl;
    } else {
      opts._tileIndexes.push_back(gridIndex);
    }

    pos += gridPattern.matchedLength();
  }

  pos = 0;
  while ((pos = coordPattern.indexIn(grids, pos)) != -1) {
    int gridIndex = coordPattern.cap(1).toInt();

    if (gridIndex < 1) {
      std::cerr << "WARNING: Ignoring tile with index of " << gridIndex << std::endl;
    } else {
      GridCoordinate coord;
      coord._x = coordPattern.cap(2).toInt();
      coord._y = coordPattern.cap(3).toInt();
      
      if (opts._gridPositions.insert(grid_coord_pair_t(gridIndex, coord)).second) {
        opts._tileIndexes.push_back(gridIndex);
        
        if (coord._x < gridBounds._gridXMin) {
          gridBounds._gridXMin = coord._x;
        }
    
        if (coord._x > gridBounds._gridXMax) {
          gridBounds._gridXMax = coord._x;
        }
    
        if (coord._y < gridBounds._gridYMin) {
          gridBounds._gridYMin = coord._y;
        }
    
        if (coord._y > gridBounds._gridYMax) {
          gridBounds._gridYMax = coord._y;
        }
      }
    }
    
    pos += coordPattern.matchedLength();
  }

  // Sort and throw away duplicates.
  std::sort(opts._tileIndexes.begin(), opts._tileIndexes.end());
  opts._tileIndexes.resize(std::unique(opts._tileIndexes.begin(), opts._tileIndexes.end()) - opts._tileIndexes.begin());

  if (opts._tileIndexes.size() == 0) {
    std::cerr << "ERROR: You must specify at least one tile to dowload." << std::endl;
    return 1;
  }

  return 0;
}

///////////////////////////////////////////////////////////////////////////////
/// \brief Parses and validates the traverse type.
///////////////////////////////////////////////////////////////////////////////
int parseTraverse(struct arg_str* parent) {
  if (strcmp(parent->sval[0], "parallel") == 0) {
    opts._zigzag = false;
  } else if (strcmp(parent->sval[0], "zigzag") == 0) {
    opts._zigzag = true;
  } else {
    return 1;
  }

  return 0;
}

///////////////////////////////////////////////////////////////////////////////
/// \brief Validates the baud rate.
///////////////////////////////////////////////////////////////////////////////
int parseBaud(struct arg_int* parent) {
  switch (parent->ival[0]) {
    case 600:
    case 1200:
    case 2400:
    case 9600:
      return 0;
    default:
      return 1;
  }
}

///////////////////////////////////////////////////////////////////////////////
/// \brief Parses and validates the traverse type.
///////////////////////////////////////////////////////////////////////////////
int parseCorner(struct arg_str* parent) {
  if (strcmp(parent->sval[0], "tl") == 0) {
    opts._corner = TOP_LEFT;
  } else if (strcmp(parent->sval[0], "tr") == 0) {
    opts._corner = TOP_RIGHT;
  } else if (strcmp(parent->sval[0], "bl") == 0) {
    opts._corner = BOTTOM_LEFT;
  } else if (strcmp(parent->sval[0], "br") == 0) {
    opts._corner = BOTTOM_RIGHT;
  } else {
    return 1;
  }

  return 0;
}

///////////////////////////////////////////////////////////////////////////////
/// \brief Parses and validates the traverse type.
///////////////////////////////////////////////////////////////////////////////
int parseDirection(struct arg_str* parent) {
  if (strcmp(parent->sval[0], "cw") == 0) {
    opts._clockwise = true;
  } else if (strcmp(parent->sval[0], "acw") == 0) {
    opts._clockwise = false;
  } else {
    return 1;
  }

  return 0;
}

///////////////////////////////////////////////////////////////////////////////
/// \brief Parses options from the command line.
///////////////////////////////////////////////////////////////////////////////
bool parseOptions(int argc, char* argv[]) {
  struct arg_lit* help     = arg_lit0(NULL, "help",           "print this help and exit");
  struct arg_str* comment  = arg_str0("m",  "comment",  NULL, "comment to be applied to each output file's metadata");

  // Surveying options.
  struct arg_str* tiles    = arg_str1("g",  "grids",    NULL, "tile(s) to download from the device");
  struct arg_int* width    = arg_int1("w",  "width",    NULL, "tile width in metres");
  struct arg_int* length   = arg_int1("l",  "length",   NULL, "tile length in metres");
  struct arg_int* numsamp  = arg_int0("s",  "numsamp",  NULL, "samples per metre (defaults to 1)");
  struct arg_int* numtrav  = arg_int0("t",  "numtrav",  NULL, "traverses per metre (defaults to 1)");
  struct arg_str* traverse = arg_str0("T",  "traverse", NULL, "traverse type [parallel|zigzag] (defaults to zigzag)");
  struct arg_str* corner   = arg_str0("c",  "corner",   NULL, "corner from which each grid was started [tl|tr|bl|br] (defaults to top left)");
  struct arg_str* direction= arg_str0("d",  "direction",NULL, "direction of the first traverse [cw|acw] (defaults to clockwise)");

  // Output options.
  struct arg_dbl* nodata   = arg_dbl0("n",  "nodata",   NULL, "dummy value to be written to the output file (defaults to -4095.0)");
  struct arg_lit* tiff     = arg_lit0(NULL, "tiff",           "output a tiff file");
  struct arg_lit* xyz      = arg_lit0(NULL, "xyz",            "output an XYZ CSV file");
  struct arg_lit* raw      = arg_lit0(NULL, "raw",            "output data as it comes from the RM15-D");

  // Georeferencing options.
  struct arg_dbl* originx  = arg_dbl0(NULL, "orgx",     NULL, "the x-coordinate of the tile origin in the projected coordinate system");
  struct arg_dbl* originy  = arg_dbl0(NULL, "orgy",     NULL, "the y-coordinate of the tile origin in the projected coordinate system");
  struct arg_dbl* axisx    = arg_dbl0(NULL, "basex",    NULL, "the x-coordinate of the x-axis in the projected coordinate system");
  struct arg_dbl* axisy    = arg_dbl0(NULL, "basey",    NULL, "the y-coordinate of the x-axis in the projected coordinate system");
  struct arg_str* srs      = arg_str0(NULL, "srs",      NULL, "The projection in either EPSG code, proj.4 format, or WKT");
  struct arg_int* forgx    = arg_int0(NULL, "forgx",    NULL, "The false origin x offset. Specified in tiles");
  struct arg_int* forgy    = arg_int0(NULL, "forgy",    NULL, "The false origin y offset. Specified in tiles");

  // Resistivity meter options.
  struct arg_int* freq     = arg_int0("f",  "freq",     NULL, "voltage frequency (Hz) (defaults to 137)");
  struct arg_dbl* probe    = arg_dbl0(NULL, "probe",    NULL, "mobile probe separation (m) (defaults to 0.5)");

  // RS232 device options.
  struct arg_int* baud     = arg_int0("b",  "baud",     NULL, "baud rate [600|1200|2400|9600]  (defaults to 9600)");
  struct arg_str* port     = arg_str1("p",  "port",     NULL, "serial port name");

  struct arg_file* outfile = arg_file1(NULL,NULL,       NULL, "output file");

  struct arg_end* end      = arg_end(20);

  void* argtable[] = {help, comment,
    tiles, width, length, numsamp, numtrav, traverse, corner, direction, 
    nodata, tiff, xyz, raw, 
    originx, originy, axisx, axisy, srs, forgx, forgy,
    freq, probe,
    baud, port,
    outfile,
    end};

  if (arg_nullcheck(argtable) != 0) {
    std::cerr << "FATAL: Insufficient memory to allocate parameters." << std::endl;
    return false;
  }

  // Defaults
  numsamp->ival[0] = 1;
  numtrav->ival[0] = 1;
  traverse->sval[0] = "zigzag";
  traverse->hdr.checkfn = (arg_checkfn*)parseTraverse;
  corner->sval[0] = "tl";
  corner->hdr.checkfn = (arg_checkfn*)parseCorner;
  direction->sval[0] = "cw";
  direction->hdr.checkfn = (arg_checkfn*)parseDirection;
  baud->ival[0] = 9600;
  baud->hdr.checkfn = (arg_checkfn*)parseBaud;
  nodata->dval[0] = -4095.0;
  freq->ival[0] = 137;
  probe->dval[0] = 0.5;
  tiles->hdr.checkfn = (arg_checkfn*)parseGrids;

  int nerrors = arg_parse(argc, argv, argtable);

  if (help->count > 0) {
    std::cout << "rm15dump version 1.0 Copyright (c) 2009, John Donovan - mersey.coder@gmail.com" << std::endl;
    std::cout << "Usage: rm15dump [options] <file>" << std::endl;
    std::cout << "Dumps data from a Geoscan Research RM15-D resistance meter." << std::endl << std::endl;
    arg_print_glossary(stdout, argtable, "  %-25s %s\n");
    std::cout << std::endl;
    std::cout << "Example usage:" << std::endl;
    std::cout << "rm15dump -g\"1,2\" -w 30 -l 30 -p \"/dev/ttyS0\" --xyz --tiff outfile" << std::endl << std::endl;
    std::cout << "Downloads the first two tiles of data from the RM15-D from the serial port /dev/ttyS0, " << std::endl;
    std::cout << "sets the width and length of the tile to 30 metres, and outputs an XYZCSV file (outfile.xyz) " << std::endl;
    std::cout << "for importing into a third-party application, and a GeoTiff (outfile.tif) file for " << std::endl;
    std::cout << "importing into a GIS application." << std::endl;

    return false;
  }

  if (nerrors > 0) {
    arg_print_errors(stdout, end, "rm15dump");
    std::cout << "Try 'rm15dump --help' for more information." << std::endl;
    return false;
  }

  if (tiff->count == 0 && raw->count == 0 && xyz->count == 0) {
    std::cerr << "ERROR: No output type selected. Use one or more of: --tiff --xyz --raw" << std::endl;
    return false;
  }

  opts._xsize = width->ival[0];
  opts._ysize = length->ival[0];
  opts._xsamp = numsamp->ival[0];
  opts._ysamp = numtrav->ival[0];
  opts._prefix = outfile->filename[0];
  opts._baud = baud->ival[0];
  opts._port = port->sval[0];
  opts._tiff = tiff->count > 0;
  opts._raw = raw->count > 0;
  opts._xyz = xyz->count > 0;
  opts._frequency = freq->ival[0];
  opts._nodata = nodata->dval[0];
  opts._probeSeparation = probe->dval[0];
  opts._comment = comment->sval[0];
  opts._originX = originx->dval[0];
  opts._originY = originy->dval[0];
  opts._axisX = axisx->dval[0];
  opts._axisY = axisy->dval[0];
  opts._coordinateRefSystem = srs->sval[0];
  opts._falseOriginX = forgx->ival[0];
  opts._falseOriginY = forgy->ival[0];

  return true;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
void fixTiles(std::vector<Sample*>& sampleVec) {
  const size_t samplesPerTraverse = opts._xsize * opts._xsamp;
  const size_t traversesPerTile = opts._ysize * opts._ysamp;
  const size_t samplesPerTile = samplesPerTraverse * traversesPerTile;
  const size_t numTiles = sampleVec.size() / samplesPerTile;

  for (size_t tile = 0; tile < numTiles; ++tile) {
    std::vector<Sample*>::iterator tileBeginIt = sampleVec.begin() + (samplesPerTile * tile);
    std::vector<Sample*>::iterator tileEndIt = sampleVec.begin() + (samplesPerTile * (tile + 1));

    // Reverse alternate traverses.
    if (opts._zigzag) {
      for (size_t i = samplesPerTraverse; i < samplesPerTile; i += samplesPerTraverse * 2) {
        std::reverse(tileBeginIt + i, tileBeginIt + i + samplesPerTraverse);
      }
    }

    // Deal with starting corner and direction.
    if (opts._clockwise) {
      switch (opts._corner) {
        case TOP_LEFT:
          // Do nothing.
          break;

        case TOP_RIGHT:
          mirrorTile(tileBeginIt, tileEndIt, samplesPerTraverse);
          flipTile(tileBeginIt, samplesPerTraverse, traversesPerTile);
          break;

        case BOTTOM_LEFT:
          flipTile(tileBeginIt, samplesPerTraverse, traversesPerTile);
          mirrorTile(tileBeginIt, tileEndIt, samplesPerTraverse);
          break;

        case BOTTOM_RIGHT:
          std::reverse(tileBeginIt, tileEndIt); // Rotate 180 degrees.
          break;
      }
    } else {
      switch (opts._corner) {
        case TOP_LEFT:
          flipTile(tileBeginIt, samplesPerTraverse, traversesPerTile);
          break;

        case TOP_RIGHT:
          mirrorTile(tileBeginIt, tileEndIt, samplesPerTraverse);
          break;

        case BOTTOM_LEFT:
          mirrorTile(tileBeginIt, tileEndIt, samplesPerTraverse);
          std::reverse(tileBeginIt, tileEndIt); // Rotate 180 degrees.
          break;

        case BOTTOM_RIGHT:
          flipTile(tileBeginIt, samplesPerTraverse, traversesPerTile);
          std::reverse(tileBeginIt, tileEndIt); // Rotate 180 degrees.
          break;
      }
    }
  }
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
void mirrorTile(std::vector<Sample*>::iterator begin, std::vector<Sample*>::iterator end, size_t samplesPerTraverse) {
  std::vector<Sample*>::iterator tileIt = begin;

  while (tileIt < end) {
    std::reverse(tileIt, tileIt + samplesPerTraverse);
    tileIt += samplesPerTraverse;
  }
}

///////////////////////////////////////////////////////////////////////////////
// Flip along TL/BR diagonal.
///////////////////////////////////////////////////////////////////////////////
void flipTile(std::vector<Sample*>::iterator begin, size_t samplesPerTraverse, size_t traversesPerTile) {
  std::vector<Sample*> buffer;
  const size_t samplesPerTile = samplesPerTraverse * traversesPerTile;

  for (size_t i = 0; i < samplesPerTile; ++i) {
    buffer.push_back(begin[((i % traversesPerTile) * samplesPerTraverse) + (i / traversesPerTile)]);
  }

  for (size_t i = 0; i < samplesPerTile; ++i) {
    begin[i] = buffer[i];
  }
}
